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ABSTRACT 

Emitter localization is a very important communications tool that will be 
extremely valuable to a multitude of different military as well as civilian applications. 
In many parts of the world, GSM is the preferred method of modulation used in 
mobile phone traffic. This thesis addresses the time difference of arrival estimation 
applied to GSM type signals using wavelet-based techniques. Signals are generated 
using the Hewlett-Packard Advanced Design System software and processed using 
algorithms based on Matlab. The results of this thesis prove that we can improve upon 
the localization of a GSM emitter through the use of wavelet-based denoising 
techniques. 
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I. INTRODUCTION 



A. BACKGROUND 

The ability to obtain localization information from a cellular phone transmission 
is increasingly becoming a much sought-after facility. Thanks to some improved 
techniques in cellular communications, this ability is very much realizable. As a matter 
of fact, recently there has been legislation which requires the development of this 
technology. When a 911 call is made from a wire line phone, the dispatcher 
automatically can provide the location of where the call was made. Currently, there is no 
such capability for calls made from cellular phones. In October 1999, President Clinton 
signed a law which makes 911 the official nationwide emergency number for both 
regular and cellular phones [Appendix A]. This legislation followed a move by the FCC 
which requires, by the end of 2001, that cellular 911 calls automatically provide the 
caller’s location [Appendix A]. 

There are countless ways in which the development of this technology will be 
beneficial to many different public and private enterprises. If this technology is utilized 
to its full potential, it will completely change the way the entire wireless industry does 
business as well provide an invaluable resource for many military applications. Listed 
here are just a few of the benefits: 

■ Location Sensitive Billing - Network operators will be able apply charges 
depending on the location of the mobile. In this way, different rates can 
be charged for calls from the home or office. This would allow network 
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operators without copper-based public switched telephone network 
(PSTN) to offer competitive rates for calls from the home or office.[l] 

■ Intelligent Transport Systems - Traffic congestion can be determined by 
the number of cellular calls coming from a particular stretch of highway. 
Vehicles with installed maps can provide directions and alternate routes to 
travelers to avoid congestion and traffic mishaps. Also companies could 
track the location of all the vehicles in its fleet and more efficiently 
dispatch vehicles according to the traffic situation and distance to the final 
destination. 

■ Enhanced Network Performance - By monitoring the frequency of calls 
from various locations over time, cellular networks can be better 
organized to provide the optimum service for a particular coverage area. 
This includes the ability to accurately monitor the movement of the mobile 
telephone to enable the network to make better decisions on when to 
initiate call handoff to a different cell. 

■ Electronic Surveillance - Law enforcement authorities can use position 
information to track the location of a mobile station which could serve a 
multitude of purposes. [2] 

■ Automobile Security - Location information can be gained from a 
strategically placed mobile transmitter in the event of a car theft. This 
would aid tremendously to the recovery of stolen vehicles and greater rate 
of apprehension of car thieves. 
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■ Unauthorized Cellular Usage Detection - Cellular carriers have long been 
plagued by the unauthorized access to personal communications systems . 
These illegal activities would diminish drastically with the advent of 
cellular phone localization. 

B. OBJECTIVE 

There are several methods available that would allow us to implement a position 
location capability into cellular communications systems. For reasons that will become 
more apparent in the following chapters, the Time Difference of Arrival (TDOA) method 
appears to be the most logical choice. In this thesis, we will discuss the TDOA 
estimation and its application to Global System for Mobile (GSM) signals. Specifically, 
we will investigate several methods to reduce the mean squared error of the TDOA 
estimate, thereby allowing a position localization with a corresponding lower error 
probablity. 

GSM was chosen to provide the backbone for this research because of its 
worldwide acceptance as the Personal Communication Service (PCS) network of choice 
and its inherent positioning localization advantages. We will discuss this system in great 
detail in order to provide an understanding of its advantages. We will also demonstrate 
our ability to simulate GSM signals using the Hewlett-Packard Advanced Design System 
(HP-ADS) program. Finally, these generated GSM signals will be evaluated using 
Matlab-based programs to provide an indication of how to improve upon the localization 
of a GSM emitter. 



3 



c. 



RELATED WORK 



The initial research on this project was begun in 1998 by Unal Aktas, LTJG, 
Turkish Navy. His work focused on the use of the wavelet transform to increase the 
accuracy of the TDOA estimation. He investigated several different denoising techniques 
and presented a comparison of the techniques. These included wavelet denoising based 
on Donoho’s method, wavelet denoising using the Wo-So-Ching threshold, wavelet 
denoising using hyperbolic shrinkage, wavelet denoising using median filtering, 
modified approximate maximum-likelihood delay estimation, denoising based on the 
fourth order moment, and a time-varying approximate maximum-likelihood technique 
The results of his work showed that the fourth order moment, and the modified 
approximate maximum-likelihood techniques were superior to the others. For that 
reason, this thesis only makes use of those three techniques. Modified versions of the 
Matlab code for wavelet-based denoising developed by Aktas, are used in this thesis. 
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II. GSM AND THE TIME DIFFERENCE OF ARRIVAL 



A. REVIEW OF GSM 

Global System for Mobile Communications (GSM) is the European standard for 
digital mobile telephones. It was developed in 1990 and first implemented in 1991 to 
replace first generation European cellular systems, which were generally incompatible 
with each other. It has gained worldwide acceptance as the first universal digital cellular 
system with modem network features extended to each mobile user, and is a strong 
contender for Personal Communication Service (PCS) above 1800 MHz throughout the 
world. [3] 

1. GSM Channel Structure 

One characteristic of GSM that makes it so desirable is that it uses both Time 
Division Multiple Access (TDMA) and Frequency Division Multiple Access (FDMA) to 
transmit and receive information. TDMA allows for multiple conversations to take place 
simultaneously at the same frequency (channel) by rationing different time slots to each 
user. GSM uses two 25 MHz bands of the radio frequency spectrum, one for uplink 
(mobile station to base station) and the other for downlink (base station to mobile). 
FDMA divides each of these blocks into 125 channels, each with a bandwidth of 200 
kHz. TDMA allows for 8 different time slots , as shown in Figure 2.1 from Ref [1], 
with each time slot being 576.92 microseconds in duration. The grouping of these 8 time 
slots is called a frame. This combined with the mechanics of FDMA provides for a total 
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of 1 000 actual speech or data channels. Currently, there is a half-rate speech coder in the 
process of development that will double the number of channels to 2000 [4], 
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Figure 2.1: GSM Frame Structure 

Of the total number of available channels, a certain number of them are dedicated 
as control channels or logical channels. Each logical channel has a specific purpose 
ranging from transmitting system parameters, to coordinating channel usage between 
base station and mobile. The logical channels are mapped onto dedicated time slots of 
particular frames. This is executed using a burst structure, the most commonly used 
being a normal burst. Contained in the middle of this burst is a 26-bit training sequence, 
chosen for its correlation properties. This received burst is correlated with a local copy of 
the training sequence to allow for the estimation of the impulse response of the radio 
channel, which aids in the demodulation of the bits in the burst. Another important 
characteristic of training sequences in the different burst structures of GSM is that they 
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lend themselves to the measurement of time which can be used to aid in the gathering of 
time-based positioning information [1]. 

2. GMSK 

GSM uses 0.3 Gaussian Minimum Shift Keying (GMSK) as its method of 
modulation. This facilitates the use of narrow bandwidth and coherent detection 
capability. In GMSK, the rectangular pulses are passed through a Gaussian filter prior to 
their passing through a modulator. Within this filter, baseband Gaussian pulse shaping 
smoothes the phase trajectory of the MSK signal, and thus stabilizes the instantaneous 
frequency variations over time. The result of this is significantly reduced sidelobe levels 
in the transmitted spectrum. The value 0.3 represents the 3 dB-bandwidth-bit duration 
product (BT) or normalized pre-Gaussian bandwidth, which corresponds to a filter 
bandwidth of 81.25 kHz for an aggregate data rate of 270.8 Kbps. As shown in Figure 
2.2 from Ref. [4], sidelobe levels fall off dramatically with decreasing BT product. 
However, reducing the BT product increases the irreducible error rate which is a product 
of intersymbol interference (ISI) created in the Gaussian filtering process. GMSK 
sacrifices the irreducible error rate for extremely good spectral efficiency and constant 
envelope properties. This is a good tradeoff as long as the irreducible error rate is less 
than that produced by the mobile channel. 

All base stations transmit a GSM synchronization burst on the synchronization 
channel. Contained in the burst is another useful 64-bit training sequence. The 
autocorrelation function for this sequence reveals a main correlation peak with a width of 
4 bits. Compared to signals used in other time-based positioning systems, this is 
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relatively wide. The GMSK modulation technique is responsible for this increased width. 
This is significant from a time-based positioning standpoint because it reduces the timing 
resolution and hence the maximum achievable position accuracy. 




NORMALIZE 0 FREQUENCY : (f-f c )T 

Figure 2.2: Power Spectral Density of a GMSK Signal 
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3. 



Interference 



In the mobile radio path, the GSM signal will experience some ISI. The mobile 
radio channel is generally hostile in nature and can affect the signal to such an extent that 
the system performance is seriously degraded [5]. Inherent in any communication system 
is the phenomenon known as multipath. This refers to the situation where a signal 
propagating from a transmitter to a receiver can travel via several different paths 
including line of sight and one or more reflected paths. The reflected paths cause delay, 
attenuation in the amplitude, and phase shifts relative to the direct path. In some 
environments, multipath components may be delayed by up to 30 microseconds. These 
delayed components cause ISI in the received signal. [6] An example of this interference 
is illustrated in Figure 2.3 from Ref. [1]. As shown, it is a distorted correlation peak 
which can cause errors in any time measurements. The degree of distortion depends on 
the relative amplitude, delay, and phase of the received signals. 

The conventional GSM receiver contains a signal demodulator which equalizes 
the received signal to improve system performance in the presence of multipath by 
combining the information received from the different arrivals. This process is called 
adaptive equalization. The equalizer does not reject multipath signals, but combines 
them. However, for positioning reasons, the receiver needs to be able to determine the 
first arrival corresponding to the line of sight, and thus requires a multipath rejection 
algorithm [1]. 
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Figure 2.3: Comparison of Ideal and Multipath Distorted Correlation Functions for GSM 
Synchronization Burst (64-Bit Training Sequence) 

4. GSM Coverage Area Layout 

The area that a GSM network covers is divided into a number of cells, each 
served by its own base station (BS), also known as a base transceiver station. The 
number of cells served by each BS, called a cluster, is determined by the amount of 
system capacity (total available channels) required for a particular area. A larger cluster 
size is an indication of a larger capacity system. One tradeoff associated with the larger 
capacity system is that a more complex system is required to control the increased 
demand for frequency reuse. Also, more cells per cluster translates into smaller 
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individual cell size and thus, more co-channel interference for which to account, as well 
as a greater frequency of handoffs between cells. The typical cluster size is 4, 7, or 12. 

Each cell is distinguished by a unique cell identifier and is allocated one or more 
uplink/downlink frequency pairs. Multiple BS’s may be grouped together under the 
control of a base station controller (BSC). In turn, several BSC’s are usually controlled 
by a mobile service center (MSC), which handles tasks such as call routing and serves as 
the interface between the mobile network and the fixed telephone network. 

5. GSM Performance Enhancing Techniques 

In addition to increasing the number of cells per cluster, GSM uses a number of 
other methods to increase capacity. These are important because they provide favorable 
implications for positioning. One such technique is the use of sectored cells. Sectoring 
effectively breaks the cell down into groups of frequencies (channels) that are only used 
within a particular region. This is best exemplified as slicing a pie into equal segments. 
Not only does this technique offer additional information for use by a positioning system, 
but it also has the effect of reducing co-channel interference. This reduced co-channel 
interference is due to the fact that a given cell receives interference and transmits with 
only a fraction of the available co-channel cells. 

Another performance enhancing technique employed by GSM is slow frequency 
hopping; so termed because the hopping rate is low compared with the symbol or bit rate. 
The mobile radio channel is a frequency-selective fading channel, which means the 
propagation conditions are different for each individual radio frequency. Slow-frequency 
hopping is used instead of fast-frequency hopping because, in order to perform well, a 
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frequency synthesizer must be able to change its frequency and settle quietly on a new 
one within a fraction of one time slot (576.92 microseconds). Frequency hopping 
provides frequency diversity to overcome Rayleigh fading due to the multipath 
propagation environment. Frequency hopping also provides interference diversity. The 
amount of interference varies with any given channel within any given cell. A receiver 
set to a channel with strong interference will suffer excessive errors over long strings of 
bursts. Frequency hopping prevents the receiver from spending successive bursts on the 
same high-interference channel. These beneficial effects of frequency hopping 
contribute to a reduction in the minimum signal-to-noise ratio (SNR) required for good 
communications. 

B. POSITIONING TECHNIQUES 

It is important to remember that in order to comply with the FCC ruling, any 
positioning solution must work with existing phones. For this reason, most cellular 
geolocation proposals are multilateral, where the estimate of the mobile’s position is 
formed by a network rather than by the mobile itself. Hence, the unilateral proposal of 
adding a global positioning system (GPS) capability to mobile phones is not a universal 
solution since mobile telephone operators would have to replace or retrofit every mobile 
telephone. Even if it were feasible to carry out this retrofit, there are other factors that 
would make this a futile attempt. First of all, the size and cost of mobile transmitters 
would be considerably increased because of the incorporated GPS receiver. Moreover, a 
clear view of the sky is necessary to receive usable GPS data from at least three satellites, 
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which makes the system useless in areas surrounded by buildings, mountains, and other 
obstructions — common working environments for cellular communications. [7] 

There are three ways that position localization information can be determined 
using the signaling aspects inherent in the GSM specification. One is through the Angle 
of Arrival (AO A) method, which uses sector information. Second, is through the use of 
propagation time using timing advances (TA). Third, is the Time Difference of Arrival 
(TDOA) method, which offers several advantages over the previous methods. In addition 
to providing some immunity against timing errors when the source of major reflections is 
near the mobile, the TDOA method is less expensive to put in place than the AOA 
method. Also, the AOA and TOA methods are hampered by requirements for line of 
sight (LOS) signal components, whereas the TDOA method may work accurately without 
a LOS component [8]. Each method will be briefly described before taking an in-depth 
look at the TDOA 

1. Angle of Arrival (AOA) 

If two or more BS platforms are used, the location of the desired mobile in two 
dimensions can be determined from the intersection of two or more lines of bearing 
(LOB), with each LOB being formed by a radial from a base station to the mobile. 
Figure 2.4 illustrates the method known as triangulation using three different platforms. 
Here, the LOB’s from the mobile to two adjacent platforms form what is called a 
measured angle, <[>. After two measured angles are calculated, trigonometry or analytic 
geometry may then be used to deduce the location of the MS. 
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Figure 2.4: Triangulation 

Another variation of the AOA method relies upon the sectoring of cells for all of 
its positioning information. Cell sectoring involves the use of highly directive antennas 
to effectively divide the cell into multiple sections. The angle of arrival can be 
determined at the base station by electronically steering the main lobe of an adaptive 
phased array antenna in the direction of the arriving mobile signal. In this case, a single 
platform may be sufficient for AOA positioning, although typically, two closely spaced 
antenna arrays are used to dither about the exact direction of the peak incoming energy to 
provide a higher resolution measurement. [9] 

AOA measurements might provide an excellent solution for wireless transmitter 
localization were it not for the fact that this method requires signals coming from the 
mobile to the base station be from the line of sight (LOS) direction only. This is seldom 
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the case in cellular systems given that the operating environment is often heavily 
cluttered with rough terrain, tall buildings and other obstructions. 

2. Time of Arrival (TOA) 

Since GSM is a TDMA system, its successful operation depends on the ability of 
all signals to arrive at the BS at the appropriate time. And since the signals arriving at the 
BS’s originate from different distances, the time at which the signals are sent must be 
varied. This is accomplished by having each BS send a timing advance (TA) to each MS 
connected to it. The TA is the amount by which the mobile station (MS) must advance 
the timing of its transmission to ensure that it arrives in the correct time slot. Obviously, 
the TA is a measure of the propagation time, which is proportional to the distance from 
the BS to the MS. Hence, the location of the MS can be constrained to a circular locus 
centered on the BS. If the MS could somehow be forced to hand over to two more BS’s, 
it would be possible to implement a positioning scheme based on the intersection of the 
three circles. 

One problem with the TOA technique is that artificially forced handoffs can be 
made to less optimal BS’s. This degrades call quality as well as reduces system capacity. 
Also, this method requires that all transmitters and receivers in the system have precisely 
synchronized clocks, since just 1 microsecond of timing error could result in a 300 meter 
position location error [9]. The use of the TA, which is in essence a timestamp, is a 
burden that most would rather not have to deal with. Furthermore, the employment of it 
presents another problem. Under GSM specifications, the TA is reported in units of a bit 
period, which equates to a locus accuracy of 554 meters. Even this considerable amount 
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of area will increase when multipath degradations and dilution of precision (the amount 
by which errors are degraded by geometry) are taken into account. 

3. Time Difference of Arrival 

The TDOA is basically an improvement upon the TOA method. The idea behind 

TDOA is to determine the relative position of the mobile transmitter by examining the 

/ 

difference in time at which the signal arrives at multiple base station receivers, rather 
than the absolute arrival time. Therefore, each TDOA measurement determines that the 
transmitter must lie on a hyperboloid with a constant range difference between the two 
receivers. The equation for this range difference is given by 

R,J = -i) ! +K - r) 1 + (Z, -z)’ -fix y -x) 1 7(rj -y) ! TtZj -z) ! (2-1) 

where the coordinates (X, Y., Z.) and (X.,Y., Z.) represent the position-fixed receivers i 
and j, and x, y, and z represent the unknown coordinate of the target transmitter. [9] 

The location of the mobile can be estimated in two dimensions from the 
intersection of two or more hyperboloids generated using TDOA measurements from 
three or more base stations. This can be extended to three dimensions if four or more 
TDOA measurements are used. 

As in the TOA method, the TDOA requires that all fixed transmitters and 
receivers in the system have precisely synchronized clocks. This, however, corresponds 
to the timing standards already in place at cellular base station sites. Unlike the TOA 
measurements, TDOA measurements do not require a timestamp of the transmitted 
signal. 
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III. POSITION LOCALIZATION USING TDOA ESTIMATION 



The process of localizing the position of a cellular transmission occurs in two 
stages. First, an accurate estimate of the TDOA measured values must be determined 
from noisy signals that are received from the mobile. This involves some type of 
denoising technique, the most popular being wavelet denoising, which is discussed 
thoroughly in [10]. Generally speaking, this involves taking the wavelet transform of the 
received signal, thresholding the wavelet coefficients, and performing the inverse wavelet 
transform on the modified coefficients. Second, a location of position estimate is 
determined from the TDOA calculations. Both phases are very complex and subject to 
minor inconsistencies that can propagate into large, unacceptable errors. Thus, accurate 
and efficient algorithms are required for both stages of processing. 

A. TDOA COMPUTATIONS 

In order to be able to evaluate the hyperbolic range equations, Eq. (2.1), one must 
first estimate the range differences Ry, or equivalently, the TDOA t, - tj. The most 
widely accepted method for obtaining these values is the generalized cross-correlation 
method. If the transmitted signal is s(t), then the signal that arrives at one receiver, x(t), 
will be delayed and corrupted by noise and is given by the equation 

x(t) = as(t- dj) + nj(t), 0 < a < 1 (3.1) 
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where dj is the delay time from the mobile to receiver one, nj(t) is noise, and a is the gain. 
Similarly, the signal that arrives at another receiver, y(t), is given by the equation 



y(t) = Ps(t - dj) + nj(t), 0 < p < 1 



(3-2) 



where dj is the delay time from the mobile to receiver two, nj(t) is the noise, and p is the 
gain. The cross-correlation function between the two received signals is derived by 
integrating the lag product of the two signals for a sufficiently long time period, T. This 
function is given by 



This approach imposes the requirement that the receiving base stations share a precise 
time reference and reference signals, but does not impose any requirement on the signal 
transmitted by the mobile. We can improve the SNR of the TDOA estimation by 
increasing the integration interval, T. The maximum likelihood estimate of the TDOA is 
obtained from the computed cross-correlation function. This estimate is equal to the 
value of r that maximizes the cross-correlation function, Eq. (3.3). 

Another method of estimating the cross correlation function, which yields the 
same results, is to compute the estimated cross-spectral density function. An inverse 
Fourier Transform of the cross-spectral density yields the estimated cross-correlation 
function. [9] 




0 



(3.3) 
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Prior to the cross-correlation function computation, some frequency domain 
processing may be used to improve the ability to distinguish the desired signal from 
arriving multipath components. This benefit is a result of the conduciveness of frequency 
domain processing to filtering techniques which are effective in reducing the effects of 
noise and interference on the TDOA estimate. 

B. OBTAINING A POSITON LOCATION ESTIMATE 

After an accurate estimate of the TDOA measurement is obtained, it is then 
possible to determine the corresponding solution to the hyperbolic equation. This process 
is initiated by substituting the range difference estimate from Eq. (2.1) into the hyperbolic 
equation, Eq. (3.3), then solving for the Cartesian coordinates of the mobile. The 
solution obtained will not be a trivial one, since the range estimates will most likely be 
noisy or inconsistent. The solution to the hyperbolic equation is also complicated by the 
fact that the fixed location receivers are arranged in a non-uniform fashion. This may be 
accounted for by linearizing the equations through Taylor-series expansion in which only 
the first two terms are retained. However, this can result in significant position location 
errors due to geometric dilution of precision (GDOP) effects. GDOP occurs when a 
relatively small range of errors results in a large position location error due to a very large 
difference in distance between the mobile unit location on the hyperbola and the locations 
of the two receivers being used. 
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There have been several proposed solutions to reduce the above inaccuracies. 
One being an exact solution for the case where the number of equations equals the 
number of unknown coordinates. However, this solution does not account for redundant 
information. Other proposals include spherical interpolation, a divide-and-conquer-based 
solution, and more recently, a proposal, by Y. T. Chan, is a closed form solution which is 
valid for both close and distant sources [9]. Chan’s method has been shown to be the 
superior of those mentioned by virtue of a lower computational complexity and higher 
noise threshold tolerance. In fact, Chan’s method guarantees optimum performance 
around small TDOA noise region since it has been shown to achieve the Cramer-Rao 
lower bound in this region. [12] 

The essence of Chan’s method is to introduce an intermediate variable to the 
original TDOA equations, then transform the nonlinear equations relating TDOA 
estimates and the source position into a set of equations which are linear in the unknown 
parameters and the intermediate variable. A least-squares computation produces an 
initial solution. By exploiting the known relation between the intermediate variable and 
the position coordinates, a second weighted least squares computation gives the final 
solution for the position coordinates. [13] 

C. DETERMINING THE ACCURACY OF POSTION LOCATION 

MEASUREMENT 

In an effort to measure the performance of TDOA position location systems, there 
are several grading criterion at one’s disposal, all of which focus upon the accuracy of the 
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position measurement. The most obvious and most widely used measure of accuracy is 
the root mean square (rms) accuracy. This is defined as the rms deviation of the 
measured position ( x,y ) about the true position (x,y), and is given by the following 



equation: 



= ^Ei(i-xy-K,y-y?) 



(3.4) 



where E is the expectation. The rms accuracy is a function of the accuracy of the raw 
locus measurements and the geometry of the base stations relative to the mobile. There is 
a limit on the achievable accuracy imposed by errors in the locus measurements. This 
accuracy can be further degraded by the physical geometry of the link between base 
stations to the mobile, referred to earlier as the GDOP. The GDOP for a 2-dimensional 
hyperbolic position location system is defined in [9] by: 



respectively, and cr R 2 is the mean square TDOA ranging error. Ideally, it would be 
prudent to minimize GDOP. However, to do this would be infringent upon a very 
important principle in the design of mobile communication systems, where the goal is to 
optimize the availability of service to mobile telephone users. Thus, we can minimize the 
GDOP, but we will degrade the availability of service in the process. Therefore, the 
focus of the accuracy determination must be placed primarily on error in the locus 
measurement. [9] 




(3.5) 



2 2 , 

where cr x and <j y are the mean square position location errors in the x and y directions. 
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IV. METHODS TO IMPROVE TDOA ESTIMATION 



Any effort to improve the TDOA estimation must focus primarily on reducing the 
effective noise at the receivers. Earlier research presented seven different denoising 
schemes designed toward this end: 1) wavelet denoising based on Donoho’s method, 2) 
wavelet denoising using the Wo-So-Ching threshold, 3) wavelet denoising using 
hyperbolic shrinkage, 4) wavelet denoising using median filtering, 5) modified 
approximate maximum-likelihood delay estimation, 6) denoising based on the fourth 
order moment, and 7) a time-varying approximate maximum-likelihood technique. The 
results of simulation of each of these methods is shown in Figure 4.1. [10] 




SNR 

Figure 4.1: Comparison of Denoising Methods for the GSM Signal 
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It was shown that method 1 failed at low SNR’s. Experimental results from [10] 
showed that method 2, while an improvement over method 1 , looses its advantage with 
increasing carrier frequency. Additional improvements are made using methods 3 and 4, 
but as can be seen in Figure 4.1, the best results (lowest mean-square error) are obtained 
using methods 5 through 7. We will take a closer look at these last three methods. 

A. FOURTH-ORDER MOMENT WAVELET DENOISING 

Inevitably, any transmitted signal will acquire some type of additive noise before 
reaching the receiver. However, it is possible to estimate the frequencies of the corrupted 
complex sinusoidal signal. A Gaussian signal being completely characterized by its 
second order statistics and the odd order moments being equal to zero for a symmetric 
probability density function, the separation of the signal and the noise requires at least the 
use of the fourth-order moments [11]. 

To define the fourth-order moment, we first model the received signal as: 

z(n) = *(«) + /(«) (4.1) 

where x(«) is a complex zero-mean, non-Gaussian, fourth order stationary signal and 
/(«) is the noise, which is a complex zero mean-Gaussian signal independent of x(n) . 
The fourth-order moment is then [11]: 

M 4 z U -Uk- i, l-i) = E(z , ,Zj,z k \z’) (4.2) 
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It was shown in earlier research that the fourth order moment of a detail function 



which contains the signal should be greater than 3cr nJ 4 , where a nd 4 is the noise power at 
subband d r [10] 

By using this property, the wavelet coefficients that represent noise can be 
eliminated while those having a signal dependency are retained. After modifying the 
detail functions, the denoised signal is obtained by performing an inverse wavelet 
transform using the modified coefficients. 

B. MODIFIED APPROXIMATE MAXIMUM LIKELIHOOD DELAY 
ESTIMATION 

Critical to the task of source localization is the time delay estimation between 
signals received at two spatially separated sensors in the presence of noise. If we let s(n) 
represent the source signal, ni(n) and n 2 (n) represent the additive noises at the respective 
sensors, and D is the difference in arrival times at the two receivers, then the receiver 
outputs, ri(n) and r 2 (n), are estimated by: 

rj(n) = ais(n) + nj(n), n = 0,1,...,T-1 (4.3) 

r 2 (n) = ct 2 s(n - D) + n 2 (n), 0 < aj, a .2 < 1 (4.4) 

where T is the number of samples collected at each channel. 

In the approximate maximum likelihood delay estimation, after wavelet 
decomposition and prior to cross-correlation, both the channel outputs are optimally 
weighted at different frequency bands with the use of an orthogonal wavelet transform. 
This weighting is done to attenuate the noise spectral components. The scaled subband 
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components are then combined using inverse wavelet transform to construct the Modified 
AML (MAML) prefiltered signal. The orthogonal wavelet transform is attractive because 
in addition to being computationally efficient, it does not suffer from the performance 
degradation due to errors that may occur in the estimation of D used in Eq. (4.4), since in 
practice, a highly accurate estimate of D may be difficult to achieve. [12] 

The final MAML delay estimate is calculated from the location of the peak of the 
cross-correlation function of the two denoised signals. 

C. TIME-VARYING MODIFIED APPROXIMATE MAXIMUM 

LIKELIHOOD DELAY ESTIMATION 

The time-varying MAML method, as is obvious from the name, works in much 
the same way as the MAML, except that it additionally accounts for the time-varying 
characteristics of the signal. Due to the time-varying nature of some signals, the 
subbands may not contain actual signal components throughout the entire segment. 
Therefore, the MAML method may loose some of its effectiveness when applied to 
segments in the presence of transitory signals. The time-varying MAML method, in 
essence, allows for the MAML method to be applied twice to each segment, thus 
approximating a time-varying gain. 
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V. SIGNAL SIMULATION AND TEST RESULTS 



A. SIGNAL GENERATION 

All signals used in this project for test purposes were generated using the Hewlett- 
Packard Advanced Design System (HP-ADS). This provided an excellent method of 
obtaining signals that are essentially the same as would be encountered in an actual 
cellular communication receiver. 

1. Advanced Design System 

The HP-ADS program is an invaluable tool for the engineers representing many 
different design aspects such as communications, digital signal processing, electronic 
circuits, mechanical circuits, and many others. Of course, for this project, we were 
particularly interested in the communications package. This package allows the user to 
custom design any type of communications system, run simulations of the design, and 
extract data collected from the simulation. Figure 5.1 shows the system used for the 
signal simulation. A more detailed view of each of the major components of this system 
can be found in Appendix [B]. Appendix [B] also provides detailed instructions for using 
the HP-ADS program to design and simulate a communications system. 

2. GSM Signal Generation 

One parameter that must be specified for the GSM signal generation is the sample 
time. The specifications of the GSM signal were presented in Chapter II. The HP-ADS 
program specifies a symbol period of 3.7037 microseconds. The filter bandwidth is 1.2 
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MHz. This allows us to sample at a minimum frequency of 2.4 MHz without violating 
the Nyquist criterion. For test purposes, we use three different sampling frequencies. As 
shown in Table 5.1, these are: 10.8 MHz, which corresponds to a sample time of 92.592 
microseconds; 5.4 MHz, which corresponds to a sample time of 185.185 microseconds; 
and 2.7 MHz, which corresponds to a sample time of 370.370 microseconds. 
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Table 5.1: GSM Signal Parameter Combinations 

The I and Q channel outputs of each of these GSM signals are shown in Figure 
5.2. One fact that is immediately evident from Figure 5.2 is that as the sampling interval 
increases, so does the variation in the GSM signals. The reason for this is obvious if we 
refer to Table 5.1. We can see that as the sample interval increases, we can expect to 
contain more symbols within a given number of samples. For test purposes, we use a 
constant number of samples (600) for all simulations. We also realize a better correlation 
characteristic with increasing sample time as shown in Figure 5.3. Better correlation is 
demonstrated by a sharper main lobe with lower sidelobes. Figure 5.3 (d) shows the 
correlation of the Matlab generated GSM signal presented by Aktas [10]. This figure 
represents the best correlation. Plots of the power spectral densities of these signals can 
be found in Appendix [C]. 
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(b) I channel; Ts = 92.592 ns 



(a) Q charnel; Ts = 92502 ns 




Figure 5.2: HP-ADS Produced GSM Signals 
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(a) Correlation of signals x and y, ADS 92.59 nsec 



(b) Comeiation of agrais x and y, ACS 185. 185 nsec 
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(c) Correlation of signals x and y; ADS 370.37nsec (d) Correlation of signals x and y; MA7LAB 





Figure 5.3: Auto-Correlation Results of Test Signals 



B. SIMULATIONS 



In Chapter IV, we discussed three methods to improve the TDOA estimate. In 
this section, we will use these methods to test our HP -ADS generated signals and see if 
we are able to extract useful TDOA data. Bear in mind that our criterion for 
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improvement in TDOA estimation is a lower mean-squared error. Each method was tried 
using four different realizations for each of the three sampling times, for a total of twelve 
different realizations per method used. However, we will only present one set from each 
sample time here, as each set follows a similar trend. The remainder of the plots can be 
found in Appendix [D]. To simplify explanation, we shall label the data sets as follows: 

a) HP-ADS generated signal with a sample interval of 92.59 nsec. 

b) HP-ADS generated signal with a sample interval of 1 85.1 85 nsec. 

c) HP-ADS generated signal with a sample interval of 370.370 nsec. 

d) Matlab generated signal. 

We will use the Matlab generated signal as our benchmark for comparison. For each 
method we will plot the mean-squared error against SNR’s in the range of -6 dB to 20 
dB. All Matlab codes used in these simulations are provided in Appendix [E]. For 
clarity in the following discussion, we define the terms average error and percentage 
improvement. Average error is calculated by adding all non-zero error values for a 
particular realization, and dividing that sum by the total number of the non-zero elements 
in that realization. Percentage improvement, as shown in the equation below, is 
calculated by subtracting the improved value from the original value and dividing the 
difference by the original value then multiplying by 1 00. 

original value - improved value 



original value 



-xl00% = percent improvement 
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mean square error m ean square error 



1 



Fourth Order Moment 



The results using this method are shown in Figure 5.4. There are two obvious 
trends that are noticed in this figure. We can easily see that the mean-squared error 
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Figure 5.4: Results of Varying Sampling Interval for Fourth-Order Moment Technique 
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improves (decreases) as the sampling interval increases. This observation is presented 
numerically in Table 5.2. From this table, we can clearly see a decrease in error values as 
we move from left to right across each row. 

In the Matlab plot (Figure 5.4 (d)), or our benchmark plot, notice that the mean- 
squared error is zero at 6 dB and higher SNR for the Fourth-Order Moment denoised 
curve. Also notice that, as expected, the denoised curve is lower than the non-denoised 
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Table 5.2: Mean-Squared Error Values at Varying SNR for the Fourth Order Moment 
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curve at all points. This shows that the denoising is effective in producing an improved 
TDOA estimate. 

Using the HP- ADS data with a sampling interval of 92.59 nsec, case(a), the mean- 
square error at -3 dB SNR is 26.3 when fourth-order denoising is employed. This is a 
3 1 .4% improvement from the non-denoised curve. The average percentage improvement 
over all SNR’s between -3 and 20 is 45.6%. 

Using HP-ADS data with a sampling interval of 185.185 nsec, we notice an 
average 52.2% increase with fourth-order moment denoising. The mean-squared error is 
5.68 at -3 dB, which is a 78.4% improvement over the -3 dB value of case (a). 

Using HP-ADS data with a sampling interval of 370.37 nsec, case (c), we obtain 
even better results. The mean-squared error at -3 dB is now 1.96, which is a 65.5% 
improvement over case (b). The average improvement realized by using fourth-order 
moment denoising instead of the non-denoised method in this case is 69.0%. We now 
see a definite trend that as our sampling interval increases, denoising increasingly 
improves the mean squared error values. 

In the benchmark case (d), there is a maximum mean-squared error of 0.25 at -3 
dB using fourth-order denoising. The average percentage improvement of fourth-order 
moment denoising over non-denoising is 72.5%. We also notice that the mean-squared 
error is zero for all SNR above 6 dB. 

The improvements we have discussed in this section can be related to the 
correlation function plots presented in Section 5. A. There, we saw that the correlation 
function plot improved with increasing sample time. Here, we state that a higher degree 
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mean square error mean square error 



of correlation of the signal components reduces the probability of error in the TDOA 
calculations. Thus, we achieve lower mean square error values with increasing degree of 
correlation of signal components. 



2. Modified Approximate Maximum Likelihood (MAML) 

The results using this method are shown in Figure 5.5. We notice a similar trend 



as in the fourth order moment method. We end up with lower mean squared error values 
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Figure 5.5: Results of Varying Sample Times for MAML Technique 



36 



as the sampling interval increases. We also notice that these mean squared error values 
are lower than those obtained in the fourth order moment method, as evidenced in Table 
5.3. 



In case (a), there is an average improvement in mean squared error of 89.6% by 
employing MAML denoising. Recall, from Section 5.B.1, the fourth order moment 
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Table 5.2: Mean-Squared Error Values at Varying SNR for MAML Method 
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denoising method produced only an average 45.6% improvement for case (a). Thus, we 
have proved that the MAML denoising method is, in fact, superior to the fourth order 
moment. 

In cases (b) and (c), there is an average improvement of 89.5% and 88.8% 
respectively. These are in accord with the averages obtained in case (a). Thus, we have 
shown that through the use of the HP-ADS program, we can generate a simulated GSM 
signal and extract the necessary data to obtain a TDOA estimation. 

3. Time-Varying Approximate Maximum-Likelihood (TV-MAML) 

In Figure 4.1, we saw that the time-varying MAML (TV-MAML) method 
produced the best results at higher SNR. The data extracted from the HP-ADS program 
did not produce similar results as obtained in the fourth order moment technique and the 
MAML method. An unsuccessful run of the TV-MAML program is shown in Figure 5.6. 
We can see in this figure that the MAML plot outperforms the TV-MAML plot. 
Furthermore, the TV-MAML plot has a non-zero mean squared error even at very high 
SNR. This can be explained as follows. Recall from Section 4.C., the TV-MAML 
requires the MAML technique be applied twice to each segment of the signal in order to 
reduce the probability of capturing non-signal components along with the actual signal. 
Recall also, all of our simulations used data sets that were 600 samples long. Therefore, 
the degree to which a particular data set is transitory depends upon its sample interval. 
Lower sampling intervals, as described in Figure 5.2, produce data sets that contain less 
information. Thus, the potential effectiveness of the TV-MAML is not fully realized if 
the signal does not contain a relatively large amount of information. 
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TV -AML Using HP-ADS Data Ts = 370.37 nsec 




Figure 5.6: Results of Unsuccessful Run of TV-MAML Program Using Sample Interval 
of 370.37 nsec 

In order to evaluate the effectiveness of the TV-MAML method, we must use a 
higher sampling interval when collecting data from the signal. Thus, we found 740 
nanoseconds to be a suitable sample time. The signal is shown in Figure 5.7 along with 
the results from the TV-MAML program. Here, we see a definite improvement relative 
to Figure 5.6. The TV-MAML outperforms the MAML program at 0 dB (and lower 
SNR’s) by an average of 18.1%. This demonstrates that the TV-MAML works better in a 
higher noise environment than the MAML. 
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0.5 




0.45 




0.4 




0.35 




0.3 


- 


0.25 




0.2 




0.15 




0.1 




0.05 




0 




-10 



(b) TVAML Method; Ts = 740 nsec 
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Figure 5.7: (a) HP-ADS Produced GSM I-Channel Signal With Sample Interval of 740 
nsec, (b) Results of TV-MAML Using Input Signal of 740 nsec Sample Interval 
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VI. SUMMARY AND CONCLUSIONS 



A. SUMMARY 

The objectives of this thesis axe to simulate GSM signals and to evaluate the mean 
squared-error for TDOA estimation. The Hewlett-Packard Advanced Design System is a 
very powerful communications development tool and proved invaluable toward our 
objective of GSM signal generation. Appendix [F] contains a concise user’s guide to the 
HP-ADS program for performing the tasks conducted in this thesis. We were able to 
manipulate this system to provide signals with desired characteristics which would later 
be used in the Matlab environment to determine the effectiveness of three different 
denoising methods. The denoising methods were presented in earlier research. These are 
the fourth-order moment, modified approximate maximum-likelihood, and time-varying 
approximate maximum-likelihood methods. We used three signal sets generated by the 
HP-ADS program to test the performance of each of these methods. The results of these 
tests agreed with the earlier research. It showed that the MAML method is the best choice 
except when the signal contains a relatively large amount of information, in which case 
the TV-MAML tends to be better. 

In the process of testing the denoising techniques, we have demonstrated our 
ability to extract localization information from our generated GSM signals. This is 
manifested in the fact that we obtain relatively low mean squared error in our TDOA 
calculations. We also found that the error could be reduced further by increasing the 
interval over which samples are taken from the signal. 
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Although the Matlab programs provided the desired results, there were some 
inconsistencies to be noted. In the Matlab code for the TV-MAML, wavelet 
decomposition and re-composition is performed on the noise-infected input signal. At 
some point between the decomposition and re-composition, a dc offset is introduced as 
well as some slow frequency oscillation toward the end of the sequence. These issues 
will be addressed in subsequent work on this project. 

B. FUTURE WORK 

In addition to providing a solution to the oscillations and DC offset described in 
Section 6.A., it may be possible to modify the existing Matlab code to obtain even better 
mean squared error values. Follow on work should also re-examine the limitations of the 
time-varying approximate maximimum-likelihood method presented in Section 5.3. An 
extension to cases with different SNR’s is also needed to determine the performance in a 
fading environment. Also, the situation where there is a weak signal in one channel and a 
stronger signal in the other channel should be examined. This would provide insight into 
the question of what ranges of SNR are reliable. Finally, an investigation into the actual 
position location estimation using the solution of the hyperbolic equation presented in 
Chapter III should be conducted. 



43 



THIS PAGE INTENTIONALLY LEFT BLANK 



44 



APPENDIX A. USA TODAY ARTICLE ON CELLULAR 911 LEGISLATION 



iS Bfr Washington 



1 0/26/99- Updated 05:39 PM ET 

Clinton signs nationwide 911 law 

WASHINGTON (AP) - President Clinton signed legislation Tuesday 
making 91 1 the official emergency number nationwide - for both 
regular and cellular phones. 

The measure also calls for development of technology that can track 
mobile callers. 

People with wireless phones now will be able to speed responses to 
highway accidents, crimes and natural disasters,” Clinton said. 
"Getting rapid care to someone who is suffering from a heart attack or 
is involved in a car crash can mean the difference between life and 
death." 

While 91 1 is widely used as the emergency number for traditional 
phones, there are 20 different codes for wireless callers across the 
country. The changes are aimed at cutting response times for the 
crews who answer 98,000 emergency calls daily from cellular phone 
callers. 

"In my home state," said Sen. Conrad Bums, R-Mont, "three quarters 
of the deaths in rural areas are because the first responders couldn't 
get there in time." 

Health care professionals joined Bums at a Capitol Hill news 
conference to applaud the new law. 

"We have great emergency room personnel. We can do a lot for 
, accident victims if we can find them and get them there," said 
Barbara Foley of the Emergency Nurses Association. "That's what 
this legislation helps us do." 

Another provision of the act directed the Federal Communications 
Commission to help states develop emergency systems, including 
technology that can automatically locate cellular callers who have 
dialed 91 1 or been involved in an accident 

The FCC in September moved forward with plans to require that 
cellular 91 1 calls automatically provide a caller's location. Regulators 
want manufacturers to begin providing locator technology within two 
years. 
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Privacy advocates have raised concerns about potential abuse of the 
technology, which would take advantage of the Global Positioning 
System developed by the military. 

The law signed Tuesday called on regulators to establish "appropriate 
privacy protection for call location information," including systems 
that provide automatic notification when a vehicle is involved in an 
accident. 

It said that calls could only be tracked in nonemergency situations if 
die subscriber had provided written approval. "The customer must 
grant such authority expressly in advance of such use, disclosure or 
access." according to Senate documents detailing provisions of fhe 
legislation. 

An estimated 700 small and rural counties have no coordinated 
emergency service to call - even with traditional phones. The bill 
would encourage private 91 1 providers to move into those areas by 
granting the same liability protections to wireless operations that now 
are offered to wireline emergency sendee systems. 

Separately, the FCC took action earlier this year to increase the 
number of cellular calls to 91 1 that are successfully completed. The 
commission required that new analog cellular phones - not existing 
phones - be made with software that routes 91 1 calls to another 
carrier when a customers own service cannot complete the call. 

Calls sometimes aren't completed because a caller is in an area where 
his or her carrier does not have an antenna, because networks are 
overloaded or because buildings or geography block signals. 

Digital phones, of which 18.8 million now are in use, were not 
covered by the new FCC rules adopted in May because such phones 
are more complex than their analog counterparts and there is no easy 
fix for the problem. 



♦ Go to Washington news 
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APPENDIX B. HP-ADS GSM COMMUNICATION SYSTEM COMPONENTS 
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APPENDIX C. POWER SPECTRAL DENSITIES OF GSM TEST SIGNALS 

PSD COMPARISON FOR SAMPLE TIME = 92.592 nsec 



Power Spectral Density Data Set 1; Ts = 92.592 nsec Power Spectral Density Data Set 2; Ts = 92.592 nsec 





Power Spectral Density Data Set 3; Ts = 92.592 nsec Power Spectral Density Data Set 4; Ts = 92.592 nsec 





49 



PSD COMPARISON FOR SAMPLE TIME = 185.185 nsec 



Power Spectral Density Data Set 1; Ts = 185 nsec Power Spectral Density Data Set 2; Ts = 185 nsec 





Power Spectral Density Data Set 3; Ts = 185 nsec 



I 
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PSD COMPARISON FOR SAMPLE TIME = 370.37 nsec 



Power Spectral Density Data Set 1; Ts = 370.37 nsec Power Spectral Density Data Set 2; Ts = 370.37 nsec 





Power Spectral Density Data Set 3; Ts = 370.37 nsec Power Spectral Density Data Set 4; Ts = 370.37 nsec 
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APPENDIX D. MATLAB SIMULATIONS FOR EACH DATA SET 



FOURTH ORDER MOMENT & MAML FOR SAMPLE TIME = 92.592 nsec 



Results for DATA 92 setl 





xcorr without denoising 
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MSE MSE 



FOURTH ORDER MOMENT & MAML FOR SAMPLE TIME = 185.185 nsec 
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MSE MSE 



FOURTH ORDER MOMENT & MAML FOR SAMPLE TIME = 370.37 nsec 
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APPENDIX E. MATLAB CODES 



FOURTH-ORDER MOMENT DENOISING TECHNIQUE 



2 -********************************************************** 

o 

************ 

9 -********************************************************** 

o 

************ 

% Denoise_sta: Wavelet Denosing Based on The Fourth Order 
Moment 

o 

o 

% SYNTAX: y=Denoise_sta (xn, yn) 

o_ 

o 

% INPUT: xn = Received signal from first receiver 
% yn = Received signal from second receiver 

% 

% OUTPUT: y = Denoised signal Xn based on Yn statistics 

% 

% SUB_FUNC : None 
% Written by Spiros Mantis 

9 -********************'*'************************************* 

o 

************ 

%********************************************************** 

************ 

function y=Denoise_s ta (xn, yn) ; 

xyn=xcorr (xn, yn, ' biased ' ) ; 

[sigmas b]=max(xyn); 
rx=xcorr (xn, ' biased ' ) ; 
maxx=rx ( length (xn) ) ; 
ry=xcorr (yn, ' biased ' ) ; 
maxy=ry ( length (yn) ) ; 
sigmanl=maxx-sigmas ; 
sigman2=maxy-sigmas ; 

lamdax=3 . l*sigmanl , '2 ; 
lamday=3 . I*sigman2^2 ; 

nx=floor (log2 (length (xn) ) ) ; 
ny=floor (log2 (length (yn) ) ) ; 
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[cx lx] =wavedec (xn, nx, ' db4 ' ) ; 
[cy ly] =wavedec (yn, ny, ' db4 ' ) ; 

dxc= [ ] ; 
for i=l:nx 

d=detcoef (cx, lx, i) ; 
dl=length (d) ; 

A= (1/dl) *sum (d. A 4 ) ; 

if A<lamdax 

dc=zeros (1, dl) ; 
else 

dc=d; 

end 

dxc=[dc dxc] ; 

end 

a=appcoef (cx, lx, ' db4 ' , nx) ; 
al=length (a); 

A=(l/al) *sum(a. A 4) ; 

if AClamdax 

ac=zeros (1, al) ; 
else 

ac=a; 

end 

dxc=[ac dxc] ; 
xd=waverec (dxc, lx, ' db4 ' ) ; 



y=xd; 
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^********************************************************** 

************ 

S2~********************************************************** 

o 

************ 

% Modif ied_tez6 : This is a test program for wavelet 
denosing 

% based on the fourth order moment technique. 

% GSM signal is used. 

o 

"o 

% SYNTAX: Modif ied_tez6 

o 

*o 

% INPUT: None 

o. 

*o 

% OUTPUT: Mean square error versus SNR 
% 

% SUB_FUNC: Denoise_sta . m 
% Written by Spiros Mantis 

O'******************************************************' k * * * 

O 

★★★★★★★★★★★★ 

&********************************************************** 

o 



clear all 

%gsm_set; %Conf iguration variables created in memory. 
% these are: 

% Tb (= 3 . 692e-6) 

% BT (= 0.3) 

% OSR (= 4) 

% SEED (= 931316785) 

% INIT_L (= 260) 

%data=data_gen ( INIT_L) ; % this creates a binary data 
% [tx_burst, I, Q] =gsm_mod (Tb, OSR, BT, data, TRAINING) ; 

%s=I+j*Q; 
data_370_set4 
s=transpose (s2) ; 
sl=length (s) ; 

pow=(l/sl) *sum(abs (s) .^2) ; 



K=100 % number of realizations 

rand ( ' seed' ,40) ; 
f=150*rand (1, K) ; 
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delay=floor (f ) ; 

delay (1 : K/2 ) =- delay ( 1 : K/2 ) ; 



% 



delay is between -150 to +150 



n= [20 18 17 16 15 14 13 12 9 6 3 0 -3 -6] ; 

SNR=10 . A (n. /10) ; 

for k=l : length (SNR) 

oran=SNR ( k) 
for i=l : K 

x= [zeros (1, 200) s zeros (1, 224 )] ; 

y= [zeros (1, 200+delay (i) ) s zeros (1, 224-delay (i) )] ; 
randn ( ' state ' , 2* ( i+j ) ) ; 

noil_real=sqrt (pow/ (2*oran) ) * randn (1,1024) ; 
randn ( ' state ' , 3* (i+j ) ) ; 

noil_imag=sqrt (pow/ (2*oran) ) *randn (1, 1024 ) ; 
randn ( ' state ' , 4 * ( i+j ) ) ; 

noi2_real=sqrt (pow/ (2*oran) ) *randn (1, 1024) ; 
randn ( ' state ' , 5* (i+j ) ) ; 

noi2_imag=sqrt (pow/ (2*oran) ) *randn (1, 1024 ) ; 

noil=noil_real+j *noil_imag; 
noi2=noi2_real+j *noi2_imag; 

xn=x+noil; % x + noise 
yn=y+noi2; % y + noise 



% TDOA calculation with xcorr ( without denoising) 

xy=xcorr (xn, yn) ; % correlation of x and y with xcorr 
[al bl] =max (real (xy) ) ; 

erl (i ) =delay (i)-(bl-1024); 



% Fourth order moment denoise 
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X_real=Denoise_sta (real (xn) , real (yn) ) ; 
X_imag=Denoise_sta (imag (xn) , imag (yn) ) ; 
Y_real=Denoise_sta (real (yn) , real (xn) ) ; 
Y_imag=Denoise_sta (imag (yn) , imag(xn)) ; 

X=X_real+ j *X_imag ; 

Y=Y_real+j *Y_imag; 

XY=xcorr (X, Y) ; % Correlation of X and Y (denoised) 

[a2 b2 ] =max (real (XY) ) ; 

erlO (i) =delay (i ) - (b2-1024 ) ; 



end 

errorl0a370set4 (k) = (1 /length (erlO) ) *sum (erlO . A 2 ) ; 
errorlsta370set4 (k) = (1 /length (erl) ) *sum(erl. A 2) ; 

H10a370set4 (k, : )=erlO; 

Hsta370set4 (k, : ) =erl; 



end 



figure ( 6) 

k= [20 18 17 16 15 14 13 12 9 6 3 0 -3 -6] ; 

plot ( k, errorlsta370set 4 (1:14), 'o' ,k, errorl0a370set4 (1:14), ' 

x ' , k, errorlsta370set4 ( 1 : 14 ) , k, errorl0a370set4 (1:14) ) 

legend( 'xcorr without denoising ' , ' sta ' ) 

title ('4TH ORDER method with cross corellation after 

denoise; 100 realizations of data_370_set4 ' ) 

ylabel ( 'MSE' ) 

xlabel ( 'SNR' ) 

figure (7 ) 

plot ( 1 : 2047 , xy ) 

title (' Correlation Function of x and y signals w/noise' ) 

figure ( 8 ) 

plot (1: 2047, XY) 

title (' Correlation Function of x and y signals without 
noise ' ) 
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save errorl0a370set4 ; 
save errorlsta370set4 

save Hsta370set4; 
save Hl0a370set4; 



MODIFIED APPROXIMATE MAXIMUM LIKELIHOOD DENOISING 



TECHNIQUE 



2~*****i'***'k'k-k'k-k'k-k'k-k'k-k-k1(-k-k'k'k'k'k'k'k'k'k'k-k-k'k'k'k'k'k-k-k-k'k'k'k-k'k-k'k'k'k'k'k'Ji:-k'k'k 

O 



o 

% Denoise: Approximate maximum-likelihood delay estimation 
% via orthogonal wavelet transform. In this function, 

% we assumed that noise is Gaussian noise and it has 

a 

% flat freq response. We modify each detail function 

% by multipliying modified AML coefficient based on 

% the signal and noise powers. 

% 

% SYNTAX: y=amll (xn, yn) 

% 

% INPUT: xn = Received signal from first receiver 

% yn = Received signal from second receiver 

o. 

'o 

% 

% OUTPUT: 

% y = Denoised signal 

% 

% SUB_FUNC : None 
% Written by Spiros Mantis 

S-********************************************************** 

o 

•k-k'k-k-k'k'k-k'k'k-k-k 



o 

-k'k-k'k'k'k'k'k'k'k'k'k 



function y=denoise (xn, yn) ; 

xyn=xcorr (xn, yn, ' biased ' ) ; 
[sigmas b]=max(xyn); 
rx=xcorr (xn, ' biased ' ) ; 
maxx=rx (length (xn) ) ; 
ry=xcorr (yn, ' biased ' ) ; 
maxy=ry (length (yn) ) ; 
sigmanl=maxx-sigmas ; 
sigman2=maxy-sigmas ; 



63 



nx=floor (log2 (length (xn) ) ) ; 
ny=floor (log2 (length (yn) ) ) ; 

[cx lx] =wavedec (xn, nx, ' db4 ' ) ; 

[cy ly] =wavedec (yn, ny, ' db4 ' ) ; 

dxc= [ ] ; 
for i=l:nx 

d=detcoef (cx, lx, i) ; 
dl=length (d) ; 
sigmad= (1/dl) *sum(d. A 2) ; 
sigmasd=sigmad-sigmanl ; 
sigmaa (i) = (2 A (nx-1) ) *sigmasd; 
if sigmasd<=0 
wd=0; 
else 

wd=sigmasd/ (sigmanl*sigman2+sigmasd* (sigmanl+sigman2) ) ; 
end 

dc=wd*d; 
dxc=[dc dxc] ; 

end 

a=appcoef (cx, lx, ' db4 ' , nx) ; 
al=length (a); 
aux=sum (sigmaa) ; 
sigmasa= ( (2 A nx) ) *sigmas-aux; 
if sigmasa<=0 
wa=0 ; 
else 

wa=sigmasa/ (sigmanl*sigman2+sigmasa* (sigmanl+sigman2) ) 

end 

ac=wa*a; 
dxc=[ac dxc]; 

xd=waverec (dxc, lx, ' db4 ' ) ; 
y=xd ; 
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★★★**★****** 

O 

************ 

o_ 

o 

% Modified_tez5 : This is a test program for modified AML 
technique. In this 

% program, we used GSM signal. 

o. 

"O 

% SYNTAX: Modif ied_tez5 

O. 

"6 

% INPUT: None 

o 

*o 

% OUTPUT: Mean square error versus SNR 

o 

"o 

% SUB_FUN: Denoise. m 
% Written by Spiros Mantis 

9 -***************★***★*************★************************ 
o 

★ ★★★★★★★★★★•A- 



clear all 



gsm_set ; 
%these are 



% Configuration variables created in memory. 

% Tb (= 3 . 692e-6) 

% BT (= 0.3) 

% OSR (= 4) 

% SEED (= 931316785) 

% INIT L (= 260) 



data=data_gen ( INIT_L) ; % this creates a binary data 
[tx_burst , I, Q] =gsm_mod (Tb, OSR, BT, data, TRAINING) ; 



s=I+ j *Q; 
sl=length (s) ; 

pow= (1/sl) *sum (abs (s) . A 2 ) ; 



K=100 % number of realizations 

rand ( ' seed' ,40) ; 
f=150*rand ( 1, K) ; 
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delay=f loor ( f ) ; 

delay ( 1 : K/2) =-delay ( 1 : K/2 ) ; 



delay is between -150 to +150 



o 

"o 



n= [20 18 17 16 15 14 13 12 9 6 3 0 -3 -6] ; 

SNR=10. A (n. /10) ; 

for k=l : length (SNR) 

oran=SNR (k) 
for i=l : K 

x= [zeros (1, 200) s zeros (1, 224 )] ; 

y= [zeros (1, 200+delay (i) ) s zeros (1,224-delay (i) )] ; 
randn ( ' state ' , 2* (i+ j ) ) ; 

noil_real=sqrt (pow/ (2*oran) ) *randn (1, 1024 ) ; 
randn ( ’ state ' , 3* (i+j ) ) ; 

noil_imag=sqrt (pow/ (2*oran) ) *randn (1, 1024 ) ; 
randn ( ' state ' , 4* (i+j ) ) ; 

noi2_real=sqrt (pow/ (2*oran) ) *randn (1, 1024 ) ; 
randn ( ' state ' , 5* ( i+j ) ) ; 

noi2_imag=sqrt (pow/ (2*oran) ) *randn (1, 1024 ) ; 

noil=noil_real+ j *noil_imag; 
noi2=noi2_real+j *noi2_imag; 

xn=x+noil; % x + noise 
yn=y+noi2; % y + noise 

% TDOA calculation with xcorr ( without denoising) 

xy=xcorr (xn, yn) ; % correlation of x and y with xcorr 
% (with the presence of noise) 

[al bl] =max (real (xy) ) ; 

erl (i) =delay (i)-(bl-1024) ; 

%AML 

x_real=denoise (real (xn) , real (yn) ) ; % Denoising of the 
%received signals 
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y_real=denoise (real (yn) , real (xn) ) ; 
x_imag=denoise (imag (xn) , imag (yn) ) ; 
y_imag=denoise (imag (yn) , imag (xn) ) ; 

X=x_real+j *x_imag; 

Y=y_real+ j *y_imag; 

XY=xcorr (X, Y) ; % Correlation of X and Y (denoised) 

[a2 b2 ] =max (real (XY) ) ; 

er8 (i) =delay (i) - (b2-1024 ) ; 



end 

error 8 a ( k) = ( 1 /length (er8 ) ) *sum (er8 . * 2 ) ; 
error la ( k) = ( 1 /length (erl) ) *sum(erl. A 2) ; 

H8a (k, : ) =er8 ; 

Hla ( k, : ) =erl ; 



end 

figure ( 6) 

k= [20 18 17 16 15 14 13 12 9 6 3 0 -3 -6] ; 

plot ( k, error la (1:14),' o' ,k, err or 8 a (1:14), 'x',k, error la (1:14 
) , k, error8a (1:14) ) 

legend (' xcorr without denoising ' , ' ami ' ) 

title ( ' Denoise AML method with cross corellation after 

denoise; 100 realizations') 

ylabel ( ' MSE ' ) 

xlabel ( ' SNR' ) 

figure (7 ) 

plot (1:2047, xy) 

title (' Correlation Function of x and y signals w/noise') 

figure ( 8 ) 

plot (1:2047, XY) 

title (' Correlation Function of x and y signals without 
noise ' ) 

save error8a; 
save errorla; 

save Hla; 
save H8a; 
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TIME- VARYING MODIFIED APPROXIMATE MAXIMUM LIKELIHOOD 



DENOISING TECHNIQUE 



2 -********************************************************** 

O 

************ 

£********************************************************** 

O 

************ 

% aml2: We modified the MAML technique by dividing each 
% detail function into two segments. Different 

% coefficients for each segment are computed. 

Q. 

O 

% SYNTAX: y=aml2 (xn, yn, delay ) 

% 

% INPUT: xn = Received signal from first receiver 

% yn = Received signal from second receiver 

% delay = True TDOA between xn and yn 

% OUTPUT: y = Error between true TDOA and estimated TDOA 

o. 

"o 

% SUB_FUNC : None 
% Written by Unal Aktas 

9 -********************************************************** 

O 

■k-Jc-k-tr-k-Jc'ki h 'k -^r -k 'k 

£********************************************************** 

O 

************ 



function y=aml2 (xn, yn, delay ) ; 

xyn=xcorr (xn, yn, 'biased' ) ; 
[sigmas b]=max(xyn); 
rx=xcorr (xn, ' biased ' ) ; 
maxx=rx ( length (xn) ) ; 
ry=xcorr (yn, ' biased ' ) ; 
maxy=ry (length (yn) ) ; 
sigmanl=maxx-sigmas; 
s igman2=maxy- sigma s ; 

nx=floor (log2 (length (xn) )) ; %10 
ny=floor (log2 (length (yn) ) ) ; 
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[cx lx] =wavedec (xn, nx, ' db4 ' ) ; %wavelet decomposition of xn 
%at level 10 

[cy ly] =wavedec (yn, ny , ' db4 ' ) ; 

dxc= [ ] ; 
for i=l:nx 

d=detcoef (cx, lx, i ) ; %detail coefficients at level i 
dl=length (d) ; %7 

Nsl=128/2 A (i-1) ; % the length of the subblock 

if Nsl<=l %i>=7; l<i>10 

N£=dl ; %7 
else 

Ns=Nsl ; % . 25<Nsl>128 

end 

D=ceil (dl/Ns) ; 
if dl<Ns*D 

dm=[d zeros (1, D*Ns-dl) ] ; 

end 

for k=l : D 

p=(k-l)*Ns+l:k*Ns; 
sigmad= (1/Ns) *sum(dm(p) . A 2) ; 
sigmasd=sigmad-sigmanl ; 
if sigmasd<=0 
wd=0 ; 
else 

wd=sigmasd/ (sigmanl*sigman2+sigmasd* ( sigmanl+sigman2 ) ) ; 
end 

dc (p) =wd*dm (p) ; 

end 

dxc=[dc (l:dl) dxc] ; 

end 

a=appcoef (cx, lx, ' db4 ' , nx) ; 
al=length (a) ; 
sigmaa= (1/al) *sum(a. A 2) ; 
sigmasa=sigmaa-sigmanl / 
if sigmasa<=0 
wa=0 ; 
else 

wa=sigmasa/ (sigmanl*sigman2+sigmasa* (sigmanl+sigman2 ) ) 

end 

ac=wa*a; 
dxc=[ac dxc] ; 
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xd=waverec (dxc, lx, ' db4 ' ) ; 



dyc= [ ] ; 
for i=l:ny 

dy=detcoef (cy, ly, i) ; 
dyl=length (dy ) ; 

Nsl=128/2 A (i-1) ; % the length of the subblock 

if Nsl<=l 
Ns=dyl; 
else 

N s =N s 1 ; 

end 

D=ceil (dyl/Ns) ; 
if dyl<Ns*D 

dmy=[dy zeros ( 1 , D*Ns )] ; 

end 

for k=l:D 

p= ( k-1 ) *Ns+l : k*Ns ; 
sigmady= (1/Ns) *sum (dmy (p) . A 2) ; 
sigmasdy=sigmady-sigmanl ; 
if sigmasdy<=0; 

wdy=0; 

else 

wdy=sigmasdy/ ( sigmanl*sigman2+sigmasdy* (sigmanl+sigman2) ) ; 
end 

dcy (p) =wdy*dmy (p) ; 

end 

dyc= [ dcy ( 1 : dy 1 ) dy c ] ; 

end 

ay=appcoef (cy, ly, ' db4 ' , ny) ; 
ayl=length (ay) ; 
sigmaay= ( 1/ayl) *sum(ay. A 2) ; 
sigmasay=sigmaay-sigmanl; 
if sigmasay<=0 
way=0; 
else 

way=sigmasay/ ( sigmanl *sigman2+sigmasay* ( sigmanl+sigman2 ) ) ; 
end 

acy=way*ay; 
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dyc=[acy dye]; 
yd=waverec (dye, ly, ' db4 ' ) ; 

% [wdl wd2 wd3 wd4 wd5 wd6 wd7 wd8 wd9 wc9] 
rxyd=xcorr (xd, yd) ; 

[a5 b5] =max (rxyd) ; 

er5=delay- (b5-1024 ) ; 
y=er5; 
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$-*************************************************- ********* 
o 

O 

************ 

% tez5t : This is a test program for time varying MAML 
technique . 

% In this program, GSM signal is used 

o. 

o 

% SYNTAX: tez5t 
% 

% INPUT: None 

a 

o 

% OUTPUT: Mean square error versus SNR 

o. 

'o 

% SUB_FUNC : aml2.m 
% Written by Unal Aktas 

%********************************************************** 

o 

************ 

°~********************************************************** 

o 

************ 



clear all 



gsm_set; 
these are: 



% Configuration variables created in memory, 

% Tb (= 3 . 692e-6) 

% BT(= 0.3) 

% OSR (= 4) 

% SEED (= 931316785) 

% INIT L (= 260) 



data=data_gen ( INIT_L) ; % this creates a binary data 
[tx_burst , I , Q] =gsm_mod (Tb, OSR, BT, data, TRAINING) ; 



%s=I+j*Q; 

data_740_setl 

s=transpose (s2) 

sl=length (s2 ) ; 

pow= ( 1/sl) *sum (abs (s ) .* 2 ) ; 



K=100 % number of realizations 

rand ( ' seed' ,40) ; 
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f=150*rand ( 1 , K) ; 
delay=f loor (f ) ; 

delay ( 1 : K/2) =-delay ( 1 : K/2) ; % delay is between -150 to 

+ 150 



n= [20 18 17 16 15 14 13 12 9 6 3 0 -3 -6] ; 

SNR=10 . A (n. /10) ; 

for k=l : length (SNR) 

oran=SNR (k) 
for i=l : K 

x= [zeros (1, 259) s zeros (1, 224) ] ; 

y=[zeros (1,259+delay (i) ) s zeros ( 1 , 224-delay ( i ))] ; 
randn ( ' state ' , 2* (i+j ) ) ; 

noil_real=sqrt (pow/ (2*oran) ) *randn (1,1024) ; 
randn ( ' state ' , 3* (i+j ) ) ; 

noil_imag=sqrt (pow/ (2*oran) ) *randn (1, 1024) ; 
randn ( ' state ' , 4* (i+j ) ) ; 

noi2_real=sqrt (pow/ (2*oran) ) *randn (1, 1024) ; 
randn ( ' state ' , 5* (i+j ) ) ; 

noi2_imag=sqrt (pow/ (2*oran) ) *randn (1,1024) ; 

noil=noil_real+ j *noil_imag; 
noi2=noi2_real+ j *noi2_imag; 

xn=x+noil; % x + noise 
yn=y+noi2; % y + noise 



% TDOA calculation with xcorr ( without denoising) 

xy=xcorr (xn, yn) ; % correlation of x and y with xcorr 
[al bl] =max (real (xy) ) ; 

erl (i) =delay (i) - (bl-1024) ; 

%AML 
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el=aml2 (real (xn) , real (yn) , delay (i) ) ; 
e2=aml2 (imag (xn) , imag (yn) , delay (i) ) ; 

ela=amll (real (xn) , real (yn) , delay (i) ) ; 
e2a=amll (imag (xn) , imag (yn) , delay (i) ) ; 

er8 (i) = (el+e2) /2; 

er8a (i) = (ela+e2a) /2; 



end 

error 8 t ( k) = ( 1 /length (er8 ) ) *sum (er8 . A 2) ; 
error la ( k) = (1/length (erl) ) *sum (erl . A 2) ; 

error 8 a ( k) = ( 1 /length (er8a) ) *sum (er8a . A 2) ; 

H8t (k, : ) =er8; 

Hla (k, : ) =erl; 



end 

%load error8a 
figure (6) 

k= [20 18 17 16 15 14 13 12 9 6 3 0 -3 -6] / 

plot ( k, err or la (1:14),' o', k, error 8 t (1:14), ' x ' , k, error 8 a (1:14 
) / ' d' , . . . 

k, error la ( 1 : 14) , k, error8t (1 : 14) , k, error 8 a (1:14) ) 

legend (' xcorr without denoising ' , ' time varying ami ' , ' ami ' ) 
title ( ' TVAML Method ADS input data; Ts = 740nsec Set#l; 100 
realizations ' ) 
figure (7) 

plot (1:541, real (s2) ) 

title ('ADS signal (s) real component; 740ns') 

save error8t; 
save errorla; 

save Hla; 
save H8t; 
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APPENDIX F. GUIDE FOR USING GSM SIMULATION IN HEWLETT- 



PACKARD ADVANCED DESIGN SYSTEM 



STEP 1 : Start Program 

If there is no icon on startup window: 
Start => Programs => HP- ADS => ADS 



STEP 2 : Open GSM Example Project 
From ADS (main) window: 

File => Example Project => Com_Sys => gsm_prj => ok => wait a few seconds 
for windows to pop up. (two windows come up: GSM_SYS Schematic^ and MODEM 
Schematic: 1 . The Modem Schematic is for a more detailed simulation. I have been primarily 
working with the GSM_SYS schematic.) DO NOT MAKE ANY CHANGES TO EITHER OF 
THESE FILES! ! ! ! Any changes will be saved in the system and they will remain there for future 
users of the program. Save to your H: drive or other drive of your choice and make changes there. 

STEP 3 : Modifying the system 

The options here are unbounded. You can experiment with all the different components in the 
Component Library. But for now, we just want to know how to get data out of this system. 

This is done by placing ‘data sinks’ in the desired locations. 

Go to the common component box at the top of the schematic window and click the 
arrow. This will show all the different components available. => 

For data extraction, you want to use ‘Sinks.’ From there you can find the appropriate 
‘sink’ for your data. (Timed and Spect. Analyzer are the most commonly used)=> 

Click once on the sink of your choice => 

Move the cursor into the schematic area and place the sink in the desired position by left- 
clicking when you have it in a clear spot. => 

End the command by using the right-click button on the mouse then click ‘end 
command’ or just hit the escape key. => 

Use wire to connect it to the system, (there is a wire button on the third row of toolbar, or 
you can go to the component button at the top) => 

Notice the name of the sink, (this is listed under the device type) This is important, 
because you will need to recall this name when you want to plot or export the data. You 
can change the name by clicking once on it and then typing the new name right there in 
the same spot. 
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STEP 4: Simulations 



After you have your desired system configuration, are ready to simulate the system: 

Simulate(on top toolbar) => simulate => wait for “simulation finished” in summary window. The 
simulation is now complete and you can extract your data. 

STEP 5 : Displaying Data 

Window(on top toolbar) => New Data Display => wait for data display window 
To see an x-y plot: 

Click once on the “plot” icon => move cursor into white area of display(a box appears), 
position the box and left-click. => choose what data you want to plot. => ok => your plot 
appears in the data display window. 

To see data: 

Click once on the “123456” icon => follow the same procedure as above. 

STEP 6 : Extracting and Exporting Data 

I haven’t figured out how to export the plots but here is how to extract the data: 

In the schematic window: click on “start the instrument server” button(third button third row of 
toolbar) => Click the following: 

Read/Write: Write 

Write to: File 

File Format: Citifile => 

Use the browse button to choose a directory and file name. => 

Choose desired data from Datasets box. => 

Click ‘write to file’ button. => your data is now saved in the specified directory. 

You will have to use wordpad or notepad to view this data. Then you will have to sort through all 
the data given to get to the specific data. This may be annoying, but rest assured the data you want 
is in there somewhere!! 



STEP 7: Importing your ADS data into Matlab 

After you have located your data in wordpad you can cut and paste from there into Matlab. 

If you want to know how many data points you have it is best to use the Matlab Editor/Debugger 
(unless you have a few hours to spend counting). 

- Make the data into a vector in Matlab and use the size command to determine the 
length. 

- Using this method you can cut it down to the right size. 
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